'''
Created on Mar 27, 2011

@author: oabalbin
'''
import sys
from collections import defaultdict, deque

def list_of_names(inputfile, pfield):
    
    genlist=deque()
    for line in inputfile:
        if line.startswith('#'):
            continue 
        fields = line.strip('\n').split('\t')
        genlist.append(fields[pfield])
    
    return genlist

def dict_of_names(inputfile, pfield):
    
    genlist=defaultdict()
    for line in inputfile:
        if line.startswith('#'):
            continue 
        fields = line.strip('\n').split('\t')
        genlist[fields[pfield]]=fields[pfield]
    
    return genlist


def read_bed_basic(ifile,entrez=True):
    '''
    chr1    16785286        16786634        NM_001145277|NECAP2     0       +
    '''
    gene_list=defaultdict()
    for l in ifile:
        f = l.strip('\n').split('\t')
        name = f[3].split('|')
        if entrez:
            gene_list[name[1]] = f
        else:
            gene_list[name[0]] = f
    
    return gene_list
    
def read_gene_base(ifile, qfile, ofile, efile, entrez=True):
    '''
    '''
    pfield=0
    ifile = open(ifile)
    qfile = open(qfile)
    ofile = open(ofile,'w')
    efile = open(efile,'w')

    gdict=dict_of_names(qfile, pfield)
    
    for l in ifile:
        if l.startswith('#'):
            continue
        f = l.strip('\n').split('\t')

        if entrez:
            gene = f[10]
        else:
            gene = f[0]
            
        try:
            loc = gdict[gene]
            ofile.write(l)
        except:
            print >> sys.stderr,"Gene not found %s"%gene
            efile.write(gene+'\n')

        
    pass

    


def genes_interest_build(exome_bed_file, gene_list_file, ofile, efile, entrez=True):
    '''
    It uses a exome bed file with gene names ref seq or ENTREZ
    and a file of genes. Output bed for gene of interest
    '''
    pfield=0
    base_ifile, query_ifile = open(exome_bed_file), open(gene_list_file)
    query_gene_list = list_of_names(query_ifile, pfield) 
    ref_gene_list = read_bed_basic(base_ifile,entrez)
    ofile = open(ofile,'w')
    efile = open(efile,'w')
    
    for g in query_gene_list:
        try:
            loc = ref_gene_list[g]
            ofile.write(",".join(loc).replace(',','\t')+'\n')
        except:
            print >> sys.stderr,"Gene not found %s"%g
            efile.write(g+'\n')
        

bfile = '/exds/projects/alignment_indexes/exome_targets/agilent/SureSelect_30Mb_Kitv2/refGene.B37.exons.SS30M.pad50.bed'
#cfile = '/exds/users/oabalbin/sw/vcfCodingSnps.v1.5/geneLists/refGene.B37.txt'
qfile = '/exds/users/oabalbin/projects/exomes/informative_genes.txt'
ofile= '/exds/users/oabalbin/projects/exomes/informative_genes_annotation.bed'
efile='/exds/users/oabalbin/projects/exomes/informative_genes_not_found.txt'
genes_interest_build(bfile, qfile,ofile,efile)
#read_gene_base(cfile, qfile, ofile, efile, entrez=True)
    
    
    